library(sp)
library(tidyverse)
library(readxl)
library(rgeos)
library('raster')
library('geosphere')
library(maptools)
library(maps)
library(haven)
options(repos=c(CRAN="http://mirror1.csvd.census.gov/CRAN"))
#install.packages('rgdal')
require ('rgdal')

# Program gets buffers around each county centroid that includes neighboring county centroids at increasing distances.

# this data is each county-year, reduce to just counties: 
setwd("//projects//")
cntyestabs <- read_dta("users//########//CountyLevelWindEstabs.dta")
cntyestabs <- subset(cntyestabs, select = c(ctyfips,x,y))
cntycenters <- na.omit(cntyestabs[!duplicated(cntyestabs$ctyfips,cntyestabs$x,cntyestabs$y),])

# number of meters in each d-mile radius
maxdist_5 <- 8047 
maxdist_10 <- 16093.5
maxdist_20 <- 32187    
maxdist_40 <- 64374
maxdist_60 <- 96561
maxdist_80 <- 128748
maxdist_100 <- 160935

cnty<-na.omit(data.frame(cbind(cntycenters$x, cntycenters$y))) 
names(cnty)<-c("Lon_C","Lat_C") 

# convert counties to spatial object
# NOTICE data.frame= keeps some data (ctyfips) for each county!
c.spoint.p <- SpatialPointsDataFrame(cnty[, c("Lon_C", "Lat_C")], 
                                     data.frame(ID=seq(1:nrow(cnty)),cntycenters[,c(1)]), 
                                     proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
# transform to planar projection: 
c.spoint.mp <- spTransform(c.spoint.p,CRS( "+init=epsg:2163" )) 

# create buffers at each distance band
for (k in c(20,40,60,80,100)) {
  #cbuf.k.p <- gBuffer(c.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE)
  assign(paste0("cbuf.",k,".p"),gBuffer(c.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE))
  assign(paste0("c.",k,".p.reproj"), spTransform(get(paste0("cbuf.",k,".p")),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
  assign(paste0("c.",k,".p.reproj"), as(get(paste0("c.",k,".p.reproj")), "SpatialPolygonsDataFrame"))  
}

writeOGR(c.20.p.reproj, dsn= "/projects/programs", layer= '20mileCirclesC', driver= 'ESRI Shapefile')
writeOGR(c.40.p.reproj, dsn= "/projects/programs", layer= '40mileCirclesC', driver= 'ESRI Shapefile')
writeOGR(c.60.p.reproj, dsn= "/projects/programs", layer= '60mileCirclesC', driver= 'ESRI Shapefile')
writeOGR(c.80.p.reproj, dsn= "/projects/programs", layer= '80mileCirclesC', driver= 'ESRI Shapefile')
writeOGR(c.100.p.reproj, dsn= "/projects/programs", layer= '100mileCirclesC', driver= 'ESRI Shapefile') 

#for (i in 1:nrow(cbuf.20.p@data)) {
# plot(cbuf.20.p[i,],add=TRUE)
#}
